/* This function stores the interpolated values in "*up" and returns 1. If two x's are very close, it returns -1 and "up" is undefined. */
# include <stdio.h>
# include <stdlib.h>

extern short interpol_1dim(double *x, double *u, unsigned short nxpt, double xp, double *up)
 {
  unsigned short i, m, ns;
  double den, dif, dift, ho, hp, w;
  double *c, *d;

  c = (double *)malloc(sizeof(double)*nxpt);
  d = (double *)malloc(sizeof(double)*nxpt);

  ns = 0;
  c[0] = d[0] = u[0];
  dif = ((dif = xp - x[0]) < 0 ? -dif:dif);

  for(i = 1; i < nxpt; i++)
       {
        if((dift = ((dift = xp - x[i]) < 0 ? -dift:dift)) < dif)
             {
              ns = i;
              dif = dift;
             }
        c[i] = d[i] = u[i];
       }

  *up = u[ns];

  for(m = 1; m < nxpt; m++)
     {
      for(i = 0; i < nxpt-m; i++)
         {
          ho = x[i] - xp;
          hp = x[m+i] - xp;
          w = c[i+1] - d[i];

          if((den = ho-hp) == (double)0) return -1; /* ERROR: Two points are very close. */

          den = w/den;
          d[i] = hp*den;
          c[i] = ho*den;
         }

      *up += (2*ns < nxpt-m ? c[ns]:d[--ns]);
     }

   free(c); free(d);

   return 1;
  }
